#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <sys/time.h>
#include <math.h>
#include <string.h>
#include <unistd.h>

// Type permettant de simuler un tableau 2D
typedef int col_t[3];

// Processus parent
#define ROOT 		0

// Informations pour le tableau d'envoi/reception et pour le tableau temporaire de voisins
#define START 		0
#define LENGTH 		1
#define RANK 		2

// Informations pour le tableau temporaire des voisins
#define STARTCOL 	2
#define STARTROW 	3
#define ENDCOL 		4
#define ENDROW 		5

#define SLEEP 		5

int rank, mpi_err,numnodes,nbproc;
int* recvcnts;
int* displs;
col_t* communicationInfos;
double* subMatrix;
double* neighborSummedValues;
int subMatrixSize = 0, nbNeighbors = 0, startIndex, startCol, startRow, endCol, endRow;
MPI_Comm comm;

// Méthode permettant de permuter deux variables
void permute ( int* val1, int* val2 )
{
    int temp = *val1;
    *val1 = *val2;   
    *val2 = temp;    
}  

// Méthode permettant à chaque CPU de connaitre ses voisins et les valeurs transmises avec ces derniers
void partition ( int n, int m )
{   
    // Pour ignorer les 0 sur la bordure de la matrice
    n-=2;
    m-=2;
 
    int childInfo[nbproc][6];
    int i,k,col,row,neighbor;
    int size = n * m, position = 0;
    int remaining = size % nbproc;   //Nombre de cases de trop à répartir entre les processeurs.
    int* partitionMatrix = malloc ( m*n*sizeof(int) );
    int tmpNeighbors [nbproc][2];
    
    subMatrixSize = size / nbproc;

    // Calcul de la position et de la longueur de chaque bloc dans la matrice
    for ( k = 0; k < nbproc; k ++ )
    {
        childInfo[k][START] = position;

        childInfo[k][LENGTH] = subMatrixSize;
        // Pour ajouter un élément à ce processus si le nombre de cases n'est pas parfaitement divisible
        if ( k < remaining )
            childInfo[k][LENGTH] ++;                
  
        // Pour calculer le debut et la fin du bloc (colonnes et rangees) dans la matrice      
        if ( childInfo[k][LENGTH] != 0 )
        {
            childInfo[k][STARTCOL] = position / n;
            childInfo[k][STARTROW] = position % n;
            childInfo[k][ENDCOL] = ( position + childInfo[k][LENGTH] - 1) / n;
            childInfo[k][ENDROW] = ( position + childInfo[k][LENGTH] - 1) % n;
            
            // Remplissage du tableau de partitionnement
            for ( i = 0; i < childInfo[k][LENGTH]; i ++ )
                partitionMatrix[childInfo[k][START]+i] = k;
        }
                
        position += childInfo[k][LENGTH];
    }    

    // Initialisation de recvcnts et de displs pour le gatherv de ROOT        
    if ( rank == ROOT )
    {
        recvcnts = malloc ( nbproc*sizeof(int) );
        displs = malloc ( nbproc*sizeof(int) ); 
         
        for ( k = 0; k < nbproc; k ++ )
        {   
            displs[k] =   childInfo[k][START];
            recvcnts[k] = childInfo[k][LENGTH];
        }
    }
    
    // Récupération des informations propres au processus courant
    subMatrixSize = childInfo[rank][LENGTH];
    startIndex = childInfo[rank][START];
    startCol = childInfo[rank][STARTCOL];
    startRow = childInfo[rank][STARTROW];
    endCol = childInfo[rank][ENDCOL];
    endRow = childInfo[rank][ENDROW];
    
    // Initialisation du tableau temporaire des voisins
    for ( k = 0; k < nbproc; k ++ )
        tmpNeighbors[k][LENGTH] = 0;
       
     
    // Calcul des voisins à gauche et à droite
    position = childInfo[rank][START] - 1;
        
    for ( i = 0; i < childInfo[rank][LENGTH]; i ++ )
    {
        col = (++position) / n;
        row = position % n;

        // Pour vérifier s'il y a un voisin à gauche
        if ( col - 1 >= 0 )
            if ( partitionMatrix[(col-1)*n + row] != rank )
            {
                neighbor = partitionMatrix[(col-1)*n+row];
                if ( tmpNeighbors[neighbor][LENGTH] == 0 )
                {
                    tmpNeighbors[neighbor][START] = i;
                    nbNeighbors ++;
                }
                tmpNeighbors[neighbor][LENGTH]++;
            }
            
        // Pour vérifier s'il y a un voisin à droite
        if ( col + 1 < m )
            if ( partitionMatrix[(col+1)*n + row] != rank )
            {
                neighbor = partitionMatrix[(col+1)*n + row];
                if ( tmpNeighbors[neighbor][LENGTH] == 0 )
                {
                    tmpNeighbors[neighbor][START] = i;
                    nbNeighbors ++;
                }
                tmpNeighbors[neighbor][LENGTH]++;
            }
    }
       
    free ( partitionMatrix );
        
    // Pour vérifier s'il y a des voisins en haut et en bas
    if ( childInfo[rank][STARTROW] != 0 )
        nbNeighbors ++;
    if ( childInfo[rank][ENDROW] < n - 1 )
        nbNeighbors ++;
           
    communicationInfos = malloc(nbNeighbors*sizeof(col_t));

    k = 0;
        
    // Pour ajouter le voisin en haut s'il y en a un (send et recv)
    if ( childInfo[rank][STARTROW] != 0 )
    {
        communicationInfos[k][START] = 0;
        communicationInfos[k][LENGTH] = 1;
        communicationInfos[k][RANK] = rank - 1;
        k++;
    }
        
    // Pour ajouter le voisin en bas s'il y en a un (send et recv)
    if ( childInfo[rank][ENDROW] < n - 1 )
    {
        communicationInfos[k][START] = childInfo[rank][LENGTH]-1; 
        communicationInfos[k][LENGTH] = 1;
        communicationInfos[k][RANK] = rank + 1;
        k++;   
    }

    // Pour ajouter tous les voisins (send et recv)
    for ( i = 0; i < nbproc; i ++ )
        if ( tmpNeighbors[i][LENGTH] != 0 )
        {
            communicationInfos[k][START] = tmpNeighbors[i][START];
            communicationInfos[k][LENGTH] = tmpNeighbors[i][LENGTH];
            communicationInfos[k][RANK] = i;
            k++;
        }
}

// Méthode permettant d'afficher la matrice generee
void displayMatrix ( int n, int m, double* matrix, char matrixRotated )
{
    int i, j;
    
    // Parce qu'une rotation a ete effectuee sur la matrice
    if ( matrixRotated )
        permute ( &n, &m );
        
    for ( j = 0; j < n; j ++ )
    {
        for ( i = 0; i < m; i ++ )
        {
            if ( i == 0 || j == 0 || i == m-1 || j == n-1 )
               printf("| %-*.3f ", 6, 0.0);
            else
                if ( matrixRotated )
                    printf("| %-*.3f ", 6, matrix[(j-1)*(m-2)+(i-1)]);
                else
                    printf("| %-*.3f ", 6, matrix[(i-1)*(n-2)+(j-1)]);                
        }
        printf("|\n");
    }    
}

// Méthode permettant d'initialiser la sous-matrice d'un processeur
void initializeSubMatrix ( int n, int m )
{
    int i, j, k;
    subMatrix = malloc ( subMatrixSize * sizeof(double) ); 
    neighborSummedValues = malloc ( subMatrixSize * sizeof(double) );
    
    for ( k = 0; k < subMatrixSize; k ++ )
    {
        i = ( startIndex + k ) / (n-2) + 1;
        j = ( startIndex + k ) % (n-2) + 1;
        subMatrix[k] = i*(m-i-1)*j*(n-j-1);
        neighborSummedValues[k] = 0;
    }    
}

// Méthode effectuant le calcul de transfert de chaleur, en parallèle
void heat(double td, double h, int n)
{
    MPI_Status status;
    MPI_Request	send_request;
    int i,j;
    double hpow = pow(h, 2);	

    // Initialisation de la matrice contenant les résultats de la deuxième partie de la
    // formule discrétisée d'Euler
    for(i = 0; i < subMatrixSize; i++)
    {
        neighborSummedValues[i] = 0.0;
    }
	
    // Envoi des valeurs aux voisins
    for(i = 0; i < nbNeighbors; i++)
    {
        MPI_Isend(&subMatrix[communicationInfos[i][START]], communicationInfos[i][LENGTH], MPI_DOUBLE, communicationInfos[i][RANK], 0, comm, &send_request);		
    }

    // Réception des valeurs des voisins
    for(i = 0; i < nbNeighbors; i++)
    {
        double recvMatrix[communicationInfos[i][LENGTH]];
        MPI_Recv(recvMatrix, communicationInfos[i][LENGTH], MPI_DOUBLE, communicationInfos[i][RANK], 0, comm, &status);

        // Calcul de la somme des valeurs voisines pour chaque cellule
        // de la sous-matrice sur ce processeur
        for (j = 0; j < communicationInfos[i][LENGTH]; j++)
        {
	    neighborSummedValues[communicationInfos[i][START]+j] += recvMatrix[j];
        }
    }

    for(i = 0; i < subMatrixSize; i++)
    {
        j = (startIndex + i) % (n-2);

        // Calcul de la somme des valeurs voisines 
        // situées dans la même sous-matrice

        // Cellule de gauche
        if ( i - (n-2) >= 0 )
        {
            neighborSummedValues[i] += subMatrix[i-(n-2)];
        }
		
        //Cellule de droite
        if ( i + (n-2) < subMatrixSize )
        {
 	    neighborSummedValues[i] += subMatrix[i+(n-2)];
        }
	
        // Cellule du haut
        if ( j - 1 >= 0 && i-1 >=0 )
        {
            neighborSummedValues[i] += subMatrix[i-1];
        }
		
	// Cellule du bas
        if ( j + 1 < (n-2) && i+1 < subMatrixSize )
        {
            neighborSummedValues[i] += subMatrix[i+1];
        }
    }
    
    for(i = 0; i < subMatrixSize; i++)
    {
        // Formule discrétisée d'Euler
        subMatrix[i] = (1.0 - 4.0 * (td / hpow)) * subMatrix[i] + (td / hpow) * neighborSummedValues[i];
        neighborSummedValues[i] = 0.0;
        usleep(SLEEP);
    }
}

// Méthode récupérant les sous-matrices de chaque processeur et les assemblant afin
// de récupérer la matrice globale.
double* rebuildMatrix ( double* subMatrix, int n, int m )
{
    double* matrix;
    
    if ( rank == ROOT )
    {
        matrix = malloc ( (n-2)*(m-2)*sizeof(double) );    
    }        

    MPI_Gatherv(subMatrix, subMatrixSize, MPI_DOUBLE, matrix, recvcnts, displs, MPI_DOUBLE, ROOT, comm);
    return matrix;    
}

// Méthode créant un communicateur pour les processeurs à utiliser pour le calcul.
void createCommGroup ( )
{
    MPI_Group g1, g2;
    int ranks [nbproc];
    int i;
  
    for ( i = 0; i < nbproc; i ++ )
        ranks[i]=i;    

    MPI_Comm_group ( MPI_COMM_WORLD, &g1 );
    MPI_Group_incl ( g1, nbproc, ranks, &g2 );
    MPI_Comm_create ( MPI_COMM_WORLD, g2, &comm );
}

// Méthode vérifiant les arguments envoyés au programme
int verifyArguments(int argc, char*argv[])
{
    char* usage = "===================================================================================\n"
                  "Exemple d'utilisation : \nmpirun --hostfile hostfile -np 24 lab3 n m tp td h nbproc\n"
                  "    n      = nombre de lignes\n"
                  "    m      = nombre de colonnes\n"
                  "    np     = nombre de pas de temps\n"
                  "    td     = temps discretise\n"
                  "    h      = taille d'un cote d'une subdivision\n"
                  "    nbproc = nombre de processus a utiliser\n"; 
    
    if ( argc != 7 )
    {
        if ( rank == ROOT )
            printf ( "Erreur : Le mauvais nombre d'arguments a ete passe au programme.\n%s", usage );
        return 0;
    }
    
    if ( atoi ( argv [1] ) < 3 || atoi ( argv [2] ) < 3 )
    {
        if ( rank == ROOT )
            printf ( "Erreur : m et n doivent etre superieurs ou egaux a 3.\n%s", usage );
        return 0;
    }

    if ( atoi ( argv [3] ) < 1 )
    {
        if ( rank == ROOT )
            printf ( "Erreur : Il doit y avoir au moins un pas de temps.\n%s", usage );
        return 0;
    }

    if ( ! ( atof ( argv [4] ) > 0 ) )
    {
        if ( rank == ROOT )
            printf ( "Erreur : td doit etre superieur a 0.\n%s", usage );
        return 0;
    }

    if ( ! ( atof ( argv [5] ) > 0 ) )
    {
        if ( rank == ROOT )
            printf ( "Erreur : h doit etre superieur a 0.\n%s", usage );
        return 0;
    }

    if ( atoi ( argv [6] ) < 1 )
    {
        if ( rank == ROOT )
            printf ( "Erreur : Il doit y avoir au moins un processeur pour effectuer le calcul.\n%s", usage );
        return 0;
    }

    return 1;
}

double euler ( double current, double left, double right, double top, double down, double td, double h )
{
    double hpow = pow ( h, 2 );
    usleep ( SLEEP );

    return (1.0-4.0*(td/hpow))*current + (td/hpow)*(left+right+top+down);
}

void runSequential ( int n, int m, int np, double td, double h )
{
    double* matrix = malloc ( n*m*sizeof(double) );
    double* oldMatrix = malloc ( n*m*sizeof(double) );
    
    int i, j, k;
    
    // Initialisation
    for ( j = 0; j < n; j ++ )
        for ( i = 0; i < m; i ++ )
            if ( i == 0 || i == m - 1 || j == 0 || j == n - 1 )
                matrix [j*m + i] = 0;
            else
                matrix [j*m + i] = i*(m-i-1)*j*(n-j-1);            

    printf ( "==== Sequentiel : matrice initiale ====\n" );
    for ( j = 0; j < n; j ++ )
    {
        for ( i = 0; i < m; i ++ )
            printf("| %-*.3f ", 6, matrix [j*m + i]);
        printf ( "|\n" );
    } 

    // Euler
    for ( k = 1; k <= np; k ++ )
    {
        memcpy ( oldMatrix, matrix, n*m*sizeof(double) );        
        
        for ( j = 1; j < n - 1; j ++ )
            for ( i = 1; i < m - 1; i ++ )
                matrix [j*m + i] = euler ( oldMatrix[j*m + i], oldMatrix[j*m + i-1], oldMatrix[j*m + i+1], oldMatrix[(j-1)*m + i], oldMatrix[(j+1)*m + i], td, h);
    }
    
    printf ( "==== Sequentiel : matrice finale ====\n" );        
    for ( j = 0; j < n; j ++ )
    {
        for ( i = 0; i < m; i ++ )
            printf("| %-*.3f ", 6, matrix [j*m + i]);
        printf ( "|\n" );
    }
    
    free ( matrix );
    free ( oldMatrix );        
}

void runParallel ( int n, int m, int np, double td, double h, int nbproc )
{
    char matrixRotated = 0;
    double* matrix;
 
    // Permuter n et m afin d'avoir des colonnes moins hautes
    if ( n > m )
    {
        permute ( &n, &m );
        matrixRotated = 1;
    }                                           

    if ( rank < nbproc )
    {    
        int i;
		
        partition ( n, m );
        initializeSubMatrix ( n, m );
		
        matrix = rebuildMatrix ( subMatrix, n, m );

        if ( rank == ROOT )
        {   
            printf ( "\n==== Parallele : Matrice initiale ====\n" );
            displayMatrix (n,m,matrix,matrixRotated);
        }
		
	for(i = 1; i <= np; i++)
	{
		heat(td, h, n);
		MPI_Barrier(comm) ;
	}
		
	matrix = rebuildMatrix ( subMatrix, n, m );
		
        if ( rank == ROOT )
        {   
		printf ( "==== Parallele : Matrice finale ====\n" );
		displayMatrix (n,m,matrix,matrixRotated);
	}
    }                                                 
}

int main(int argc,char *argv[])
{
    int n, m, np;
    double td, h;
    struct timeval tp;
    double time1, time2, time3;
        
    // Initialisation
    mpi_err = MPI_Init(&argc,&argv);
    mpi_err = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    mpi_err = MPI_Comm_size(MPI_COMM_WORLD, &numnodes);    

    if ( verifyArguments ( argc, argv ) )
    {
        n  = atoi ( argv [1] );
        m  = atoi ( argv [2] );
        np = atoi ( argv [3] );
        td = atof ( argv [4] );
        h  = atof ( argv [5] );
        nbproc = atoi ( argv [6] );

        // Parce qu'il y a moins de cases que de processeurs.
        if ( (n-2)*(m-2) < nbproc )
        {
            nbproc = (n-2)*(m-2);
            if ( rank == ROOT )  
                printf("nbproc a ete change a %d du au faible nombre de cases.\n", nbproc);
        }

        // Création d'un groupe de communication personnalisé
        // afin de n'inclure que les processus nécessaires.
        createCommGroup ();
        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
        // Chronometre pour le début du traitement parallele        
        if ( rank == ROOT )
        { 
            gettimeofday (&tp, NULL);
            time1 = (double) (tp.tv_sec) + (double) (tp.tv_usec) / 1e6;
        }

        runParallel ( n, m, np, td, h, nbproc );
    
        if ( rank == ROOT )
        {
            // Chronomètre pour la fin du traitement parallele/le début du traitement sequentiel
            gettimeofday (&tp, NULL);
            time2 = (double) (tp.tv_sec) + (double) (tp.tv_usec) / 1e6;
 
            runSequential ( n, m, np, td, h );
            
            // Chronometre pour la fin du traitement sequentiel
            gettimeofday (&tp, NULL);
            time3 = (double) (tp.tv_sec) + (double) (tp.tv_usec) / 1e6;
   
            printf("Temps d'execution parallele : %f s\n", time2 - time1 );
            printf("Temps d'execution sequentiel : %f s\n", time3 - time2 );
            printf("Acceleration : %f\n", ( time3 - time2 ) / ( time2 - time1 ) );
        }
    }
    
    MPI_Barrier( MPI_COMM_WORLD ) ;

    mpi_err = MPI_Finalize();
	
    return 0;
}
